{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "85a8fc3f",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# 导言区"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f397dd62",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import root\n",
    "from CoolProp.CoolProp import PropsSI as psi\n",
    "import scipy.constants as sc\n",
    "from Functions.UnsteadyStateConduction import *\n",
    "from Functions.Self_defined import check_Fo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeb7fc36",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题04-01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cb540e45",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.         3.625      1.73515625 1.86412598 1.98521817 2.00374495\n",
      "  2.00151149 2.0001245  1.9999501 ]\n",
      " [0.         5.675      4.54546875 4.02875293 3.97865685 3.99416033\n",
      "  3.99986273 4.00026367 4.00006125]\n",
      " [0.         3.76875    4.99605469 5.06074878 5.0127267  4.99958744\n",
      "  4.99927858 4.99987183 5.00000964]]\n"
     ]
    }
   ],
   "source": [
    "N = 9\n",
    "t = np.zeros((3, N))\n",
    "# xpr1 = 8*t[0] + 2*t[1] + t[2] - 29\n",
    "# xpr2 = t[0] + 5*t[1] + 2*t[2] - 32\n",
    "# xpr3 = 2*t[0] + t[1] + 4*t[2] - 28\n",
    "for i in range(N-1):\n",
    "    t[0][i+1] = 1/8 * (29 - 2*t[1][i] - t[2][i])\n",
    "    t[1][i+1] = 1/5 * (32 - t[0][i+1] - 2*t[2][i])\n",
    "    t[2][i+1] = 1/4 * (28 - 2*t[0][i+1] - t[1][i+1])\n",
    "print(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56903b86",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题04-02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "79037b5f",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[300.         275.         258.984375   252.70385742 250.75216293\n",
      "  250.21956414 250.06397532 250.01850614 250.00536644]\n",
      " [300.         268.75       256.54296875 251.74407959 250.50418377\n",
      "  250.14840923 250.04286684 250.01241597 250.00360406]\n",
      " [200.         167.1875     154.27246094 151.26457214 150.37407279\n",
      "  150.10749204 150.03115772 150.00904981 150.00262326]\n",
      " [200.         160.546875   153.31420898 150.99210739 150.28155893\n",
      "  150.08176405 150.02378326 150.00688899 150.00199743]]\n"
     ]
    }
   ],
   "source": [
    "t_a = 100\n",
    "t_b = 500\n",
    "t_c = 100\n",
    "t_d = 100\n",
    "\n",
    "N = 9\n",
    "t = np.zeros((4, N))\n",
    "t[0][0] = t[1][0] = 300\n",
    "t[2][0] = t[3][0] = 200\n",
    "for i in range(N-1):\n",
    "    t[0][i+1] = 1/4 * (t_a + t_b + t[2][i] + t[1][i])\n",
    "    t[1][i+1] = 1/4 * (t_b + t_c + t[0][i+1] + t[2][i])\n",
    "    t[2][i+1] = 1/4 * (t_c + t_d + t[1][i+1] + t[3][i])\n",
    "    t[3][i+1] = 1/4 * (t_d + t_a + t[2][i+1] + t[0][i+1])\n",
    "print(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb722b5d",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题04-03"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "db4d9849",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "可以使用集中参数法\n",
      "tau = 329 s\n"
     ]
    }
   ],
   "source": [
    "d = 5e-2\n",
    "l = 30e-2\n",
    "t_0 = 30\n",
    "t_oo = 1200\n",
    "t = 800\n",
    "h = 140\n",
    "c = 0.48e3\n",
    "rho = 7753\n",
    "lambda_ = 33\n",
    "\n",
    "r = d / 2\n",
    "V = np.pi * r**2 * l\n",
    "A = 2 * np.pi * r**2 + np.pi * d * l\n",
    "l_c = V / A\n",
    "Bi = get_Bi(l_c, lambda_, h)\n",
    "\n",
    "if Bi < 0.05:\n",
    "    print(f'可以使用集中参数法')\n",
    "    Delta_T_ratio = (t - t_oo) / (t_0 - t_oo)\n",
    "    a = get_a(lambda_, rho, c)\n",
    "    tau_c = get_tau_c(l_c, lambda_, a, h)\n",
    "    tau = - tau_c * np.log(Delta_T_ratio)\n",
    "    print(f'tau = {tau:.0f} s')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28321d6f",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-04"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "eb68f031",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "剖面上的最大温差为：129.66 C\n"
     ]
    }
   ],
   "source": [
    "delta = 100e-3\n",
    "t_oo = 1000\n",
    "t_0 = 20\n",
    "t = 500\n",
    "\n",
    "h = 174\n",
    "lambda_ = 34.8\n",
    "a = 0.555e-5\n",
    "Bi_list = [0.1, 0.5, 1.0]\n",
    "mu_list = [0.3111, 0.6533, 0.8603]\n",
    "shape = 'P'\n",
    "\n",
    "Bi = get_Bi(delta, lambda_, h)\n",
    "mu = np.interp(Bi, Bi_list, mu_list)\n",
    "\n",
    "x = delta\n",
    "eta = x/delta\n",
    "\n",
    "theta = t - t_oo\n",
    "t_m = t_oo + theta / theta_to_theta_m_ratio(mu, eta, shape)\n",
    "delta_T_m = abs(t - t_m)\n",
    "print(f'剖面上的最大温差为：{delta_T_m:.2f} C')\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    Fo = p\n",
    "    xpr = (t - t_oo) / (t_0 - t_oo) - theta_to_theta_0_ratio(mu, eta, Fo, shape)\n",
    "    return xpr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e90f6a3a",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "03fbc4bc",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "所需的时间为：5532 s\n"
     ]
    }
   ],
   "source": [
    "d = 400e-3\n",
    "t_0 = 20\n",
    "t_oo = 900\n",
    "t = 750\n",
    "\n",
    "h = 174\n",
    "lambda_ = 34.8\n",
    "a = 0.695e-5\n",
    "shape = 'C'\n",
    "\n",
    "l_c = d / 2\n",
    "eta = 1\n",
    "Bi = h * l_c / lambda_\n",
    "mu = get_mu(Bi, shape)\n",
    "\n",
    "ratio_w_to_0 = (t - t_oo)/(t_0 - t_oo)\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    Fo = p\n",
    "    return ratio_w_to_0 - theta_to_theta_0_ratio(mu, eta, Fo, shape)\n",
    "\n",
    "\n",
    "guess_values = np.ones(1)\n",
    "Fo = root(expressions, guess_values).x[0]\n",
    "tau = Fo * l_c**2 / a\n",
    "print(f'所需的时间为：{tau:.0f} s')\n",
    "\n",
    "check_Fo(Fo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fc64485",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-06"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "26c17ab5",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 706.08 C\n"
     ]
    }
   ],
   "source": [
    "t_0 = 20\n",
    "t_w = 1450\n",
    "x = 80e-3\n",
    "tau = 2 * 3600\n",
    "a = 0.89e-6\n",
    "\n",
    "t = t_x_for_constant_t_w(x, tau, t_0, t_w, a)\n",
    "print(f't = {t:.2f} C')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdb026b4",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-07"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9682d101",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = 0.87 m\n"
     ]
    }
   ],
   "source": [
    "t_0 = 10    # 原地表温度，能否作为原地表以下的土壤的统一温度，有疑问\n",
    "t_w = -15\n",
    "tau = 45 * sc.day\n",
    "t = 0\n",
    "c = 1840\n",
    "rho = 2050\n",
    "lambda_ = 0.52\n",
    "\n",
    "a = lambda_ / (rho * c)\n",
    "\n",
    "guess_values = np.ones(1)\n",
    "x = root(lambda x: t - t_x_for_constant_t_w(x, tau, t_0, t_w, a), guess_values).x[0]\n",
    "print(f'x = {x:.2f} m')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d377bcf6",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-08"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7fe17afe",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "最低温度为：1157.66 C\n",
      "最高温度为：1197.97 C\n"
     ]
    }
   ],
   "source": [
    "shape = 'P'\n",
    "delta_1 = 0.5/2\n",
    "delta_2 = 0.7/2\n",
    "delta_3 = 1/2\n",
    "lambda_ = 40.5\n",
    "a = 0.722e-5\n",
    "t_oo = 1200\n",
    "tau = 4*3600\n",
    "t_0 = 20\n",
    "h = 348\n",
    "\n",
    "Bi_1 = h * delta_1 / lambda_\n",
    "Fo_1 = a * tau / delta_1**2\n",
    "mu_1 = get_mu(Bi_1, shape)\n",
    "ratio_w_to_0_1 = theta_to_theta_0_ratio(mu_1, 1, Fo_1, shape)\n",
    "ratio_m_to_0_1 = theta_to_theta_0_ratio(mu_1, 0, Fo_1, shape)\n",
    "\n",
    "\n",
    "Bi_2 = h * delta_2 / lambda_\n",
    "Fo_2 = a * tau / delta_2**2\n",
    "mu_2 = get_mu(Bi_2, shape)\n",
    "ratio_w_to_0_2 = theta_to_theta_0_ratio(mu_2, 1, Fo_2, shape)\n",
    "ratio_m_to_0_2 = theta_to_theta_0_ratio(mu_2, 0, Fo_2, shape)\n",
    "\n",
    "\n",
    "Bi_3 = h * delta_3 / lambda_\n",
    "Fo_3 = a * tau / delta_3**2\n",
    "mu_3 = get_mu(Bi_3, shape)\n",
    "ratio_w_to_0_3 = theta_to_theta_0_ratio(mu_3, 1, Fo_3, shape)\n",
    "ratio_m_to_0_3 = theta_to_theta_0_ratio(mu_3, 0, Fo_3, shape)\n",
    "\n",
    "\n",
    "ratio_m_to_0 = ratio_m_to_0_1 * ratio_m_to_0_2 * ratio_m_to_0_3\n",
    "t_m = t_oo + ratio_m_to_0 * (t_0 - t_oo)\n",
    "print(f'最低温度为：{t_m:.2f} C')\n",
    "\n",
    "ratio_w_to_0 = ratio_w_to_0_1 * ratio_w_to_0_2 * ratio_w_to_0_3\n",
    "t = t_oo + ratio_w_to_0 * (t_0 - t_oo)\n",
    "print(f'最高温度为：{t:.2f} C')\n",
    "\n",
    "check_Fo(Fo_1, Fo_2, Fo_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6185c849",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-09"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1c8cb647",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t_m = 1181.26 C\n"
     ]
    }
   ],
   "source": [
    "shape = ['P', 'C']\n",
    "d = 600e-3\n",
    "l = 1000e-3\n",
    "t_0 = 30\n",
    "t_oo = 1300\n",
    "tau = 4*3600\n",
    "h = 232\n",
    "lambda_ = 40.5\n",
    "a = 0.625e-5\n",
    "\n",
    "# 先考虑厚度方向，作为无限大平板进行分析\n",
    "l_c_1 = l / 2\n",
    "Bi_1 = h * l_c_1 / lambda_\n",
    "Fo_1 = a * tau / l_c_1**2\n",
    "mu_1 = get_mu(Bi_1, shape[0])\n",
    "eta_1 = 0\n",
    "ratio_m_to_0_1 = theta_to_theta_0_ratio(mu_1, eta_1, Fo_1, shape[0])\n",
    "\n",
    "# 再考虑径向，作为无限长圆柱进行分析\n",
    "l_c_2 = d / 2\n",
    "Bi_2 = h * l_c_2 / lambda_\n",
    "Fo_2 = a * tau / l_c_2**2\n",
    "mu_2 = get_mu(Bi_2, shape[1])\n",
    "eta_2 = 0\n",
    "ratio_m_to_0_2 = theta_to_theta_0_ratio(mu_2, eta_2, Fo_2, shape[1])\n",
    "\n",
    "ratio_m_to_0 = ratio_m_to_0_1 * ratio_m_to_0_2\n",
    "t_m = t_oo + ratio_m_to_0 * (t_0 - t_oo)\n",
    "print(f't_m = {t_m:.2f} C')\n",
    "\n",
    "check_Fo(Fo_1, Fo_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7bac104",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dbdb4690",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tau = 1721 s\n",
      "Q = 28067 J\n"
     ]
    }
   ],
   "source": [
    "D = 4e-2\n",
    "H = 6e-2\n",
    "shape = ['P', 'C']\n",
    "t_0 = 10\n",
    "t_oo = 180\n",
    "h = 15\n",
    "t_m = 80\n",
    "material = 'water'  # 用水的物性参数替代牛肉的物性参数\n",
    "# 以（10°C+80°C)/2=45°C来确定从开始加热到中心温度为 80c水的物理特性，\n",
    "# 以（10°C+180°C)/2=95°C来确定计算总加热量的物性。\n",
    "\n",
    "# 先求出牛肉块的质量\n",
    "T_0 = sc.convert_temperature(t_0, 'C', 'K')\n",
    "rho_0 = psi('D', 'T', T_0, 'P', sc.atm, material)\n",
    "V = np.pi * D**2 * H / 4\n",
    "m = rho_0 * V\n",
    "\n",
    "t_a = (t_0 + t_m) / 2\n",
    "t_b = (t_0 + t_oo) / 2\n",
    "\n",
    "T_a = sc.convert_temperature(t_a, 'C', 'K')\n",
    "rho_a = psi('D', 'T', T_a, 'P', sc.atm, material)\n",
    "c_a = psi('C', 'T', T_a, 'P', sc.atm, material)\n",
    "lambda_a = psi('L', 'T', T_a, 'P', sc.atm, material)\n",
    "a_a = lambda_a / (c_a * rho_a)\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    tau = p\n",
    "    # 先考虑厚度方向，作为无限大平板进行分析\n",
    "    l_c_1 = H / 2\n",
    "    Bi_1 = h * l_c_1 / lambda_a\n",
    "    Fo_1 = a_a * tau / l_c_1**2\n",
    "    mu_1 = get_mu(Bi_1, shape[0])\n",
    "    eta_1 = 0\n",
    "    ratio_m_to_0_1 = theta_to_theta_0_ratio(mu_1, eta_1, Fo_1, shape[0])\n",
    "\n",
    "    # 再考虑径向，作为无限长圆柱进行分析\n",
    "    l_c_2 = D / 2\n",
    "    Bi_2 = h * l_c_2 / lambda_a\n",
    "    Fo_2 = a_a * tau / l_c_2**2\n",
    "    mu_2 = get_mu(Bi_2, shape[1])\n",
    "    eta_2 = 0\n",
    "    ratio_m_to_0_2 = theta_to_theta_0_ratio(mu_2, eta_2, Fo_2, shape[1])\n",
    "\n",
    "    ratio_m_to_0 = ratio_m_to_0_1 * ratio_m_to_0_2\n",
    "    xpr = t_m - t_oo - ratio_m_to_0 * (t_0 - t_oo)\n",
    "    return xpr\n",
    "\n",
    "\n",
    "guess_values = np.ones(1)\n",
    "tau = root(expressions, guess_values).x[0]\n",
    "print(f'tau = {tau:.0f} s')\n",
    "\n",
    "# 先考虑厚度方向，作为无限大平板进行分析\n",
    "l_c_1 = H / 2\n",
    "Bi_1 = h * l_c_1 / lambda_a\n",
    "Fo_1 = a_a * tau / l_c_1**2\n",
    "mu_1 = get_mu(Bi_1, shape[0])\n",
    "Q_to_Q_0_1 = Q_to_Q_0_ratio(mu_1, Fo_1, shape[0])\n",
    "\n",
    "# 再考虑径向，作为无限长圆柱进行分析\n",
    "l_c_2 = D / 2\n",
    "Bi_2 = h * l_c_2 / lambda_a\n",
    "Fo_2 = a_a * tau / l_c_2**2\n",
    "mu_2 = get_mu(Bi_2, shape[1])\n",
    "Q_to_Q_0_2 = Q_to_Q_0_ratio(mu_2, Fo_2, shape[1])\n",
    "\n",
    "Q_to_Q_0 = Q_to_Q_0_1 + Q_to_Q_0_2 * (1 - Q_to_Q_0_1)\n",
    "\n",
    "T_b = sc.convert_temperature(t_b, 'C', 'K')\n",
    "rho_b = psi('D', 'T', T_b, 'P', sc.atm, material)\n",
    "c_b = psi('C', 'T', T_b, 'P', sc.atm, material)\n",
    "Q_0 = c_b * m * (t_oo - t_0)\n",
    "Q = Q_0 * Q_to_Q_0\n",
    "print(f'Q = {Q:.0f} J')\n",
    "\n",
    "check_Fo(Fo_1, Fo_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "459ed5e2",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7bbb4e80",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1) tau = 1 s\n",
      "(2) tau = 60 s\n",
      "(3) tau = 6036 s\n",
      "(4) tau = 603587 s\n"
     ]
    }
   ],
   "source": [
    "t_0 = 25\n",
    "t_w = 50\n",
    "x = np.array([0.01, 0.1, 1.0, 10])\n",
    "delta_t = 0.1\n",
    "a = 1e-5\n",
    "t = t_0 + delta_t\n",
    "\n",
    "theta_ratio = (t - t_w) / (t_0 - t_w)\n",
    "eta = root(lambda eta: theta_ratio - erf(eta), 1).x[0]\n",
    "tau = (x/2/eta)**2/a\n",
    "\n",
    "for i, tau_ in enumerate(tau):\n",
    "    print(f'({i+1}) tau = {tau_:.0f} s')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33ccb0c6",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b927416b",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "T_oo = 303\n",
    "f = 5\n",
    "v = 20\n",
    "d = 0.9e-3\n",
    "rho = 8332\n",
    "c_p = 188\n",
    "Lambda_ = 51\n",
    "\n",
    "# 涉及到对流换热公式，后续再求解。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da5e98a3",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "40096c8c",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t_0_tau = 215.44 C\n",
      "t_x_tau = 46.12 C\n"
     ]
    }
   ],
   "source": [
    "q_0 = 2e4\n",
    "t_0 = 20\n",
    "a = 1e-7\n",
    "lambda_ = 0.2\n",
    "t_w_max = 180\n",
    "tau = 30\n",
    "x = 3e-3\n",
    "\n",
    "t_0_tau = t_x_for_constant_q_0(0, tau, t_0, lambda_, a, q_0)\n",
    "t_x_tau = t_x_for_constant_q_0(x, tau, t_0, lambda_, a, q_0)\n",
    "print(f't_0_tau = {t_0_tau:.2f} C')\n",
    "print(f't_x_tau = {t_x_tau:.2f} C')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5e726d9",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 例题03-14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t_m = 98.51 C\n",
      "Q = 165211 J\n"
     ]
    }
   ],
   "source": [
    "d = 10e-2\n",
    "height = 8e-2\n",
    "t_0 = 40\n",
    "t_oo = 105\n",
    "tau = 80*60\n",
    "shape = ['P', 'C']\n",
    "\n",
    "# 物性参数可以选择水在温度为初始温度和最终温度的平均值时的物理参数\n",
    "# 相变换热系数h很大，可认为h = np.inf\n",
    "material = 'Water'\n",
    "h = np.inf\n",
    "\n",
    "# 先求出蔬菜罐头的质量\n",
    "T_0 = sc.convert_temperature(t_0, 'C', 'K')\n",
    "rho_0 = psi('D', 'T', T_0, 'P', sc.atm, material)\n",
    "V = np.pi * d**2 / 4 * height\n",
    "m = rho_0 * V\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    t = p\n",
    "    t_a = (t_0 + t)/2\n",
    "    T_a = sc.convert_temperature(t_a, 'C', 'K')\n",
    "    rho = psi('D', 'T', T_a, 'P', sc.atm, material)\n",
    "    c = psi('C', 'T', T_a, 'P', sc.atm, material)\n",
    "    lambda_ = psi('L', 'T', T_a, 'P', sc.atm, material)\n",
    "    a = lambda_ / (rho * c)\n",
    "\n",
    "    # 先考虑厚度方向，作为无限大平板进行分析\n",
    "    l_c_1 = height / 2\n",
    "    Bi_1 = h * l_c_1 / lambda_\n",
    "    Fo_1 = a * tau / l_c_1**2\n",
    "    mu_1 = get_mu(Bi_1, shape[0])\n",
    "    eta_1 = 0\n",
    "    ratio_m_to_0_1 = theta_to_theta_0_ratio(mu_1, eta_1, Fo_1, shape[0])\n",
    "    Q_to_Q_0_1 = Q_to_Q_0_ratio(mu_1, Fo_1, shape[0])\n",
    "\n",
    "    # 再考虑径向，作为无限长圆柱进行分析\n",
    "    l_c_2 = d / 2\n",
    "    Bi_2 = h * l_c_2 / lambda_\n",
    "    Fo_2 = a * tau / l_c_2**2\n",
    "    mu_2 = get_mu(Bi_2, shape[1])\n",
    "    eta_2 = 0\n",
    "    ratio_m_to_0_2 = theta_to_theta_0_ratio(mu_2, eta_2, Fo_2, shape[1])\n",
    "    Q_to_Q_0_2 = Q_to_Q_0_ratio(mu_2, Fo_2, shape[1])\n",
    "\n",
    "    ratio_m_to_0 = ratio_m_to_0_1 * ratio_m_to_0_2\n",
    "    t_m = t_oo + ratio_m_to_0 * (t_0 - t_oo)\n",
    "    xpr = t - t_m\n",
    "    return xpr\n",
    "\n",
    "\n",
    "guess_value = t_0 + 10\n",
    "t = root(expressions, guess_value).x[0]\n",
    "t_m = t\n",
    "print(f't_m = {t_m:.2f} C')\n",
    "\n",
    "t_a = (t_0 + t)/2\n",
    "T_a = sc.convert_temperature(t_a, 'C', 'K')\n",
    "rho = psi('D', 'T', T_a, 'P', sc.atm, material)\n",
    "c = psi('C', 'T', T_a, 'P', sc.atm, material)\n",
    "lambda_ = psi('L', 'T', T_a, 'P', sc.atm, material)\n",
    "a = lambda_ / (rho * c)\n",
    "\n",
    "# 先考虑厚度方向，作为无限大平板进行分析\n",
    "l_c_1 = height / 2\n",
    "Bi_1 = h * l_c_1 / lambda_\n",
    "Fo_1 = a * tau / l_c_1 ** 2\n",
    "mu_1 = get_mu(Bi_1, shape[0])\n",
    "eta_1 = 0\n",
    "ratio_m_to_0_1 = theta_to_theta_0_ratio(mu_1, eta_1, Fo_1, shape[0])\n",
    "Q_to_Q_0_1 = Q_to_Q_0_ratio(mu_1, Fo_1, shape[0])\n",
    "\n",
    "# 再考虑径向，作为无限长圆柱进行分析\n",
    "l_c_2 = d / 2\n",
    "Bi_2 = h * l_c_2 / lambda_\n",
    "Fo_2 = a * tau / l_c_2 ** 2\n",
    "mu_2 = get_mu(Bi_2, shape[1])\n",
    "eta_2 = 0\n",
    "ratio_m_to_0_2 = theta_to_theta_0_ratio(mu_2, eta_2, Fo_2, shape[1])\n",
    "Q_to_Q_0_2 = Q_to_Q_0_ratio(mu_2, Fo_2, shape[1])\n",
    "\n",
    "Q_to_Q_0 = Q_to_Q_0_1 + Q_to_Q_0_2 - Q_to_Q_0_1 * Q_to_Q_0_2\n",
    "t_b = (t_0 + t_m)/2\n",
    "T_b = sc.convert_temperature(t_b, 'C', 'K')\n",
    "c_b = psi('C', 'T', T_b, 'P', sc.atm, material)\n",
    "Q_0 = c_b * m * (t_oo - t_0)\n",
    "Q = Q_0 * Q_to_Q_0\n",
    "print(f'Q = {Q:.0f} J')\n",
    "\n",
    "check_Fo(Fo_1, Fo_2)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}